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ABSTRACT 

Aims. Solving the continuum radiative transfer equation in high opacity media requires sophisticated numerical tools. In 
order to test the reliability of such tools, we present a benchmark of radiative transfer codes in a 2D disc configuration. 
Methods. We test the accuracy of seven independently developed radiative transfer codes by comparing the temperature 
structures, spectral energy distributions, scattered light images, and linear polarisation maps that each model predicts 
for a variety of disc opacities and viewing angles. The test cases have been chosen to be numerically challenging, with 
midplane optical depths up 10®, a sharp density transition at the inner edge and complex scattering matrices. We also 
review recent progress in the implementation of the Monte Carlo method that allow an efficient solution to these kinds 
of problems and discuss the advantages and limitations of Monte Carlo codes compared to those of discrete ordinate 
codes. 

Results. For each of the test cases, the predicted results from the radiative transfer codes are within good agreement. The 
results indicate that these codes can be confidently used to interpret present and future observations of protoplanetary 
discs. 

Key words, radiative transfer — circumstellar matter — accretion discs — planetary systems: proto-planetary discs 
— methods: numerical 
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J> 1. Introduction lengths dust re-emits the absorbed radiation. How much 
• , radiation is scattered and absorbed is a function of both the 
^ . Dust represents an essential element in the energy balance geometry of the circumstellar environment and the proper- 
H ; of a variety of astrophysical objects, from the mterstellar ^jgg ^-^^ ^^^^^ ^^^.^^ the amount of absorbed radiation 
. , medmm to the atmospheres and close circumstellar envi- gg^-g ^^le temperature of the dust (and gas) and defines the 
ronments of numerous classes of object; from the lower mass amount of radiation that is re-emitted at longer, thermal 
planets and brown dwarfs, to massive stars. With the ad- wavelengths 
vent of high-angular resolution and high-contrast instru- 
ments, the basic structural properties (e.g., size, inclina- To get a reliable understanding of the structure and 
tion, and surface brightness) of the circumstellar environ- evolution of these "dusty" objects, be it the evolution of 
ments of the nearest and/or largest objects — discs and dust grain sizes, the temperature dependent chemistry, 
envelopes around young stars in nearby star- forming re- or simply the density profiles, it is highly desirable to 
gions and around more distant evolved stars — are now model not only the integrated fluxes {i.e., the spec- 
under close scrutiny. With this unprecedented wealth of tral energy distributions, hereafter SED), but also the 
high-resolution data, from optical to radio, detailed studies resolved brightness and/or polarisation profiles when 
of the dust content become possible and sophisticated ra- available. This can only be done by solving the radiative 
diative transfer (RT) codes are needed to fully exploit the transfer (hereafter RT) problem in media that can have 
data. large optical depths and/or complex geometries and 
At short wavelengths, dust grains efficiently absorb, compositions. Recent studies of circumstellar discs are 
scatter, and polarise the starlight while at longer wave- based on detailed comparisons of high-quality data sets, 
combining various kinds of observation (SED, multiple 
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codes (e.R IWood et all 120021: IWatson fc Stapelfeldtl l2004t 
Wolf et al.l l2003t iD'Alessio et al I !2006l: 'Steinacke r et al I 



2006t iDoucetet aLl l2007t iFitzgerald ct al 



Pontoppidan et al 



20071: iPinte et all l2007l . 



Glauser et al.ll2008i:lfMlnirkulam et al.ll2008D . Such studies 



will become more and more common with the advent of 
new instruments like VLT/SPHERE, Gemini/GPI, JWST, 
Herschel and ALMA, and validating the accuracy of RT 
codes is of particular importance. 

Analytical solutions do not exist for wavelength- 
dependent radiative transfer and sophisticated numerical 
methods must be used. Testing the reliability of RT com- 
putations requires in that case to compare the solutions to 
well-defined prob lems by independe nt codes. Such a work 
has been done bvllvezic et al.l (1 19971) f or a ID spherical ge- 
ometry and by iPascucci et al.l (|2004l hereafter P04) m a 
2D disc configuration. The later work compared in detail 
the calculations of five radiative codes and has been used 
as a reference to v alidate newly develop ed RT cod es (e 
iHarries et alll2004 lErcolano et'aLll2005i iPinte et al.|[200 
The test cases were however limited to relatively modest op- 
tical depths (midplane opacity r < 100 in the V band), or- 
ders of magnitude smaller than the actual optical depths of 
protoplanetary discs for which radiative transfer codes are 
generally used in the literature. High optical depths rep- 
resent challenging calculations for radiative transfer codes, 
with potential convergence issues, and additional tests are 
required to confidently trust the results of radia tive codes 
in this regime. Furthermore, calculations in lP04l were done 
assuming isotropic scattering and restricted to SEDs. With 
the advent of high-resolution observations of young stel- 
lar objects, validating the calculations of resolved surface 
brightness of these objects is now also crucially needed. 



In this paper, we extend the work of |P04| to realistic, 
very optically thick discs with anisotropic scattering. We 
perform a comparison of seven RT codes in a well defined 
2-dimensional disc configuration, with simple dust prop- 
erties. The prediction of four Monte Carlo codes (tem- 
perature structure in the disc, the emergent SED as well 
as monochromatic scattered light images, and polarisation 
maps for different disc opacities and viewing angles) are 
compared in the case of anisotropic scattering (sections [2] 
to HI). Additional comparisons with discrete ordinate codes, 
in the case of isotropic scattering and without scattering, 
are presented in section [S] 



2. Radiative transfer modelling 

2.1. The radiative transfer problem 

Solving the RT problem in dusty environments aims at 
determining the (polarised) specific intensity Ix{'r^,lt) at 
each point and direction Ti of the volume and at each 
wavelength A. This intensity is obtained by solving the sta- 
tionary transfer equation. 

In the case of randomly oriented dust particles, the ra- 
diative transfer equation can be written adopting the Stokes 
formalism: 



ds 



= -«r(7)JA(r,7?) 



Sx{r,lt',lt)Ix{-^,lt')An' (1) 



where I\i/r , ~n) — (/, Q, [/, V) is the Stokes vector, with / 
representing the total intensity, Q and U the linearly po- 
larised intensities, and V the circularly polarised intensity. 



^abs 



( r ), K3f^( r ) and ( r ) 



[t) + K5f^( r ) are 



the absorption, scattering and extinction opacities, respec- 
tively, s is the length along the direction of propagation. 
Sxi/r is the 4x4 scattering (or Mueller) matrix 

describing the changes in the Stokes vector when the light is 
scattered from the direction r?' to the direction la . B\{T) 
is the Planck function and /q is the unitary Stokes vector 
representing unpolarised emission Jq — (1,0,0,00 

Computation of the thermal emission requires to deter- 
mine the dust temperature structure T{~r^). This tempera- 
ture is determined by writing that the dust is in radiative 
equilibrium. If we assume that the dust is at the local ther- 
modynamic equilibrium and that there is no more sources 
of energy than the radiation field, the temperature is ob- 
tained by solving the implicit equation: 



^t%r)B,{T{y))dx= I Kr{r)j,d\ (2) 



^abs / 



where J\ is the mean specific intensity {i.e. the specific 
intensity averaged over all solid angles). 

The system of equations [T] and [2] completely defines the 
RT problem when the dust optical properties {Kf^^, k'^^^, 
Sx) and sources of radiation (initial conditions for equa- 
tion [1]) are given. It is important to note that opacities can 
depend on the temperature, in which case solving simulta- 
neous equations [1] and [2] requires an iterative scheme. Most 
of the time however, dust opacities do not vary much with 
temperature and can be assumed to be constant. We will 
make this assumption in the following analysis. 

Additionally, radiative transfer plays an integral role in 



the physics of the disc. It alt ers the density stru cture via 
hydrostatic equilibrium (e.g. I Walker et al .11200^ and im- 
pacts the dust content, and hence the opacity. For instance 
differential dust sublimation at the inner edge can destroy 
some of the dust grains, th e temperature of a dust grain 
depending on its size (e.g., ' Tannirkulam et al] 120071 ) and 
composition (e.g.. lWo itkc 200(|). Similarly the formation of 
ice mantles around the grains in the cold, outer regions of 
the disc also affects the grain opacities. Including any of 
these effects requires an iterative approach. In this paper, 
we restrict ourselves to the benchmark of radiative transfer 
solvers and keep the density structure and dust properties 
fixed. 



2.2. Numerical methods 

2.2.1. The Monte Carlo method 

Anisotropic scattering by dust grains precludes the use of 
direct methods to solve the continuum RT equation and 



^ because the grains are randomly oriented, the nett dust ther- 
mal emission is not polarised. 
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Fig. 1. Dust optical properties for the 1 /im silicate grains used in the calculations. Left: dust opacity, the full line 
represents the extinction opacity, the dashed line the absorption opacity and the dotted line the scattering opacity. 
Right: albedo (full line) and asymmetry parameter g = (cosO) (dashed line). For wavelengths larger than 100 /im, both 
the albedo and asymmetry parameter have values close to and the contribution from scattered light is negligible. 



Monte Carlo methods are commonly used instead. They 
solve the RT equation by stochastically propagating "pho- 
ton packets" through the dusty environment. The transport 
of packets is governed by scattering, absorption and re- 
emission events that are controlled by the optical properties 
of the medium (opacity, albedo, scattering phase function, 
etc) and by the temperature distribution. Upon leaving the 
model boundaries, "photon packets" are used to build an 
SED and/or synthetic images. 

The Monte Carlo scheme estimates physical quan- 
tities by statistical means, which potentially leads to 
noisy results when the number of packets sampling 
some regions and/or directions in the model becomes 
low. Several variance reduction techniques have been de- 
veloped to improve the sam pling of the Monte Carl o 
method: forced first scattering (Ca shwcU fc Iilverettlll959l ). 
peel-off techniques (jYusef-Zadeh et al.. 1981) . estimation 
of the mean specific intensity ( Lucvl IT999I ). immediate 
re-emission ([Biorkm an fc WoodI 20011 ). and importance 
weighting schemes (|juvcla 2005^. These various techniques 
allowed the Monte Carlo method to progressively become 
competitive against grid-based methods, and it is now more 
and more commonly used to solve the continuum RT prob- 
lem. 

Despite these techniques, Monte Carlo methods can be- 
come computationally expensive when the optical depth 
becomes very large. In our disc configuration, this leads to 
two major difhculties: 

— because gradients of opacity are oriented toward the disc 
surface, packets entering the disc tend to escape after 
a few interactions. Very few packets penetrate the cen- 
tral regions of the disc, leading to a noisy temperature 
structure close to the disc midplane. 

— for edge-on configurations, the flux in the near- and mid- 
infrared is dominated by packets originating from the 
dense central regions of the disc. These packets will be 
scattered towards the observer in the surface layers of 
the disc, where the probability of interaction is very low 
(due the very low density). Packets escaping the dense 



regions (after a large number of interactions) almost 
always go through the surface layers without interacting 
with the dust grains and large number of packets are 
required to converge SEDs or images in this regime. 

In the next paragraphs, we describe two schemes that, 
when coupled with a Monte Carlo approach, significantly 
improve the efficiency of the codes. They have helped to 
overcome the previously mentioned difficulties and to effi- 
ciently solve the test cases presented in this paper. 

2.2.2. Diffusion approximation 

In the deep regions of the disc, solving the complete RT 
equation is not necessary. The radiation field becomes 
isotropic and the source function becomes equal to the 
Planck function. The behaviour of the density of energy 
e(T*) = 4:aT{~r^)'^/c is in that case properly described by 
the diffusion theory: 



V. {D{T')Ve{'r)) 



(3) 



where the diffusion coefficient is defined as D{~t^) = 
l/2>p{jr)nR(jr) with K/j(T*) the local Rosseland opacity. 

In this case, Monte Carlo methods can be efficiently cou- 
pled to diffusion approximation methods. The model can be 
divided into two regions. The first one, that corresponds the 
surface of the disc represents all parts of the model volume 
where the optical depth in any direction is smaller than a 
given threshold. In this region, the temperature structure is 
computed with a Monte Carlo method, eventually includ- 
ing acceleration schemes such as the immediate re-emission 
concept and/or estimation of the mean specific intensity. 

In the rest of the model volume, i.e. in the central re- 
gions of the disc, the temperature structure can be solved 
using the diffusion approximation. The temperature struc- 
ture at the edge of the diffusion approximation region (ini- 
tial condition for equation [3]) is given by the solution of the 
Monte Carlo calculations. 

The optical depth threshold which defines these two re- 
gions must be set high enough to ensure that the radiation 
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field inside the diffusion approximation region is isotropic 
and dominated by the local emission. 

For an optimal efficiency, this method must avoid cal- 
culating the complete propagation of photon packets in- 
side the diffusion approximation region. This can be done 
by using various methods. The propagation of the pack- 
ets can be calculated in a faster w ay by using a modified 
rando m walk procedure ([Fleck fc'Ca nfield ,19841 : [Min et al.l 
I2OO9I) . This method combines multiple interaction steps in 
one computation, while still keeping track of the energy de- 
posited in an accurate sense. In this way, the interaction 
between the optically thick regions and the upper layers of 
the disc is properly computed. In addition, the number of 
steps taken can be easily adjusted to the local radiation field 
and density gradient, making it a highly flexible method. A 
"mirror" condition can also be used: when a packet enters 
the diffusion approximation region, it is sent back with the 
same energy and wavelength but with an opposite direction 
vector. Although not rigorously exact, this method (used by 
TORUS and MCFOST) was found to be very accurate when 
compared to the full Monte Carlo solution (see section HTTjl . 

Details and tests on the accu racy of diffusion approxi- 
mation methods are presented in lMin et al.l ()2009[ ). 

2.2.3. Ray-tracing 

When the specific intensity is known in each point of the 
model, direct ray-tracing methods using the formal solution 
of the RT equation can be used to produce observables. 

Ray-tracing has for instance been used in combination 
with Monte Carlo methods to produce SEDs and emis- 
sion maps in the infrared and millimetre regimes, where 
scattering c an be considered isotropic in some cases 
l2003; DuUe mond fc Domini^ 120041. When the scattering 
is isotropic, only the mean specific intensity and temper- 
ature structure are required to calculate the source func- 
tion. These two quan tities can b e easily estimated with a 
Monte Carlo method (|Lucvlll999f) and do not require large 
amount of memory to be stored. After an initial Monte 
Carlo run computing the total mean intensity and temper- 
ature structure in the disc, SEDs and/or maps can then be 
produced by integrating the source function on rays origi- 
nating from the observer. Which such a method, the Monte 
Carlo method is used only to estimate the specific intensity 
and not the images and/or SEDs. The resulting noise in the 
observables is much lower as it only refiects the noise in the 
mean specific intensity and no longer the noise associated 
to the production of the observables themselves. 

This method combining Monte Carlo and ray-tracing 
can be extended to any wavelength if the angular depen- 
dence of the scattering component of the source function 
is preserved in the calculations. The Monte Carlo method 
produces all the information needed to perform such calcu- 
lations, as it can give an estimate of the specific intensity, 
with its complete angular dependence, and not only of the 
mean specific intensity. However, storing the full spatial, 
angular and wavelength dependence of the radiation field 
requires large amounts of memory which is currently be- 
yond computational capacities. 

This difficulty can be overcome by invoking successive 
monochromatic Monte Carlo runs, which removes the need 
to store the wavelength dependence of the specific inten- 
sity. An initial multi-wavelength Monte Carlo run calcu- 
lates the temperature structure in the disc. The SED is 



then constructed wavelength by wavelength with succes- 
sive monochromatic Monte Carlo runs that estimate the 
specific intensity at each point of the model. 

In MCFOST, the specific intensity is then saved for a 
set of angular directions (method 1). At the end of each 
Monte Carlo run, the scattering emissivity in any direc- 
tion is calculated from the specific intensity and resolved 
maps and/or integrated fluxes for any inclinations are fl- 
nally obtained by ray-tracing. This step (monochromatic 
Monte Carlo run -I- ray-tracing run) is repeated over all 
wavelengths, without storing the specific intensity at the 
previous wavelength. 

A slightly different method is adopted in MCMax, where 
the scattering emissivity in a given set of directions is stored 
instead of the specific intensity itself (method 2) . The scat- 
tering emissivity is the product of the specific intensity by 
the scattering phase function, i.e. the last term in equa- 
tion [T] Each time a packet crosses a cell, its contributions 
to the scattering emissivity in the chosen directions are cal- 
culated by multiplying the packet energy by the local phase 
function. At the end of each monochromatic Monte Carlo 
run, maps and/or fiuxes at the chosen inclinations are pro- 
duced by integrating the source function via a ray-tracing 
method. 

Method 1 avoids the expensive calculations of the scat- 
tering emissivity each time a packet crosses a cell but re- 
quires a larger amount of memory to store the angular 
dependence of the specific intensity. If the radiation field 
is stored for a few specific wavelengths, it also allows to 
produce scattered light images and emission maps at any 
inclinations, by only running additional ray-tracing calcu- 
lations. However, with this method, the angular sampling 
of the radiation field must be performed with care, espe- 
cially when scattering is very anisotropic. This issue is not 
encountered with method 2, which has the same, almost 
perfec10, angular sampling of the radiation field as classical 
Monte Carlo methods. 

2.3. Codes description 
2.3.1. MCFOST 

MCFOST is a 3D continuum and line radiative transfe r 
code based on the Monte Carlo method ijPinte et al.ll2006l ). 
Temperature structure s are calculated usiiig the i mmediate 
re-emission concept of iBiorkman fc WoodI (|200ll ) but with 
a conti nuous depos ition of energy to estimate the mean in- 
tensity (!Lucvlll999D . The code uses a spherical or cylindrical 
grid, with an adaptive mesh refinement at the inner edge 
(based on the opacity gradient) so as to properly sample 
the inner radius of the disc. 

Several improvements have been i mplemented on to p 
of the original algorithm presented in iPinte et al.l (j2006[ ). 
In very optically thick parts of the model, the tempera- 
ture structure is calculated with a diffusion approximation 
method, using the Monte Carlo calculations as limit con- 
ditions. Equation [3] is solved as an asymptotic limit of the 
time dependent diffusion equation, via an implicit scheme 
to ensure stability and accuracy. The transition between the 
Monte Carlo and diffusion approximation regions is set to 
an optical depth of 1 000 at the wavelength where the stel- 
lar emission peaks. To avoid calculating the propagation of 

^ only limited by the numerical precision 
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photon packets inside the diffusion approximation domain, 
packets are mirrored at the boundaries of the Monte Carlo 
domain. 

The temperature structure and radiation field estimated 
by the Monte Carlo runs are used to produce images, polar- 
isation maps, and SEDs with a ray-tracing method, where 
the emerging flux is obtained by calculating the formal so- 
lution of the radiative equation along rays. 



2.3.2. MCMax 

The Monte Carlo radiative transfer code MCMax is based 
on the scheme of immed iate re-emission as proposed by 
IBiorkman fc WoodI (|200ll ). For the te mperature str ucture 
the method of continuous absorption bv lLucvl (|l999 f) is im- 
plemented. The photons are traced in 3D on a spherical 
coordinate grid while the geometry of the system is set 
to be cylindrically symmetric. For the optically thick re- 
gions a modified random walk procedure is applied in or- 
der to mak e multiple interaction steps in a single computa- 
tion (see iFleck fc Canfieidlll984l : lMin et al.ll2009} ). This has 
the advantage that the computational speed is increased 
significantly while the temperature structure is still com- 
puted with high accuracy. After the Monte Carlo procedure 
a partial diffusion approximation is used for these regions 
in the disc that received too few p hotons to determ ine a re- 
liable temperature structure fsee lMin et al.l [20091 ). All the 
observables are constructed by integrating the formal solu- 
tion using ray-tracing. In this way noise on the observables 
is reduced significantly. 

The spatial grid at the inner edge of the disc is set in 
such a way that the optical depth for both the local radia- 
tion field and the stellar radiation are sampled logarithmi- 
cally. 



2.3.3. Pinball 

Pinball is a Monte-Carlo code that calculates scattered-light 
images; it does not calculate the equilibrium temperature or 
include thermal re-emission. An earlier version and a pair 
of sim ple test cases were described by Watso n fc HennevI 
(|2001h . The current version includes polarisation. 



2.3.4. TORUS 

TORUS is a 3D continuum and line radi ative transfer 
code based on the Monte Carlo metho d (Harriet l2000t 
iHarries et al.ll200S iKurosawa et al.ll2004D . Radiative equi- 
librium is c omputed usi ng the continuous absorption algo- 
rithm from iLucvl ()1999[ ) . Calculations are performed on a 
2D, cylindrical adaptive-mesh grid. Storing the opacity in- 
formation on an adaptive mesh has particular advantages 
for the problem considered here, since it allows an adequate 
sampling of the inner edge of the disc, where the opacity 
gradient is very steep. The temperature structure in the 
central regions of the disc is computed with a diffusion ap- 
proximation method. 

F or sc attered light images, the enforced scattering con- 
cept (ICashwell fc Everett 1959j ) as well as the peel-off tech- 
nigue (jYusef-Zadeh et al.l 19841 ) are implemented to reduce 
the variance. 



The descriptions of the codes ProDiMo, RADMC and 
RADICAL, that have performed the test cases with isotropic 
scattering and without scattering, are presented in sec- 
tion O 



3. Benchmark problem 

All the codes in this paper that calculate the thermal 
equilibrium have successfully reproduced the |P04| bench- 
mark, from optically thin configurations to optical depths 
of 100 in the optical. The t est c ases presented here are 
complementary to those in |P04| and are restricted to 
optical depths higher than 1 000. The full description 
of the benchmark problem, including tabulated values 
for the disc density and the dust properties, as well as 
the results for all codes are presented on the webpage 
http : / / www . astro . ex . ac . uk/ people/cpinte/benchmark/j 
This should allow additional codes to compare their results 
with the ones presented in this paper. 

3.1. System geometry 

The model consists of a dusty disc surrounding a central 
star and located at a distance of 140 pc. 

We consider an axisymmetric flared den- 
sity structure with a Gaussian vertical profile 
p{r,z) = /9o(r) exp(— z^/2 /i(r)^). We use power-law 
distributions for the surface density S(r) — Sq (r/ro)^^-^ 
and the scale height h{r) = ho (r/ro)^'^^^ where r is the 
radial coordinate in the equatorial plane and ho — 10 AU 
is the scale height at the radius tq = 100 AU. The disc 
extends from an inner cylindrical radius rin = 0.1 AU to 
an outer limit rout — 400 AU. The edges of the disc are 
assumed to be sharp, i.e. vertical: there is nothing inside 
Tin and outside rout and the density is defined by the 
previously mentioned power-laws between them. The dust 
disc mass is the only parameter varied and takes 4 different 
values: 3 x 10-^ 3 x 10"^ 3 x 10"^ and 3 x IO^^Mq. 

This configuration was chosen because it represents a 
mo re di fficult problem to solve than the test case presented 
bv |P04 the disc extends much closer to the star, 0.1 AU 
instead of 1 AU and the radial gradient of density is much 
steeper with a slope of surface density of —1.5 instead of 
0.125, leading to much higher densities close to the inner 
edge of the disc, and hence much higher disc optical depths. 

The star is defined as a uniformly radiating blackbody 
sphere at a temperature of 4 000 K and with a radius of 2 
solar radii. 



3.2. Dust properties 

Dust grains are defined as homogeneous and spherical par- 
ticles with a singl e size of 1 fim and are compo sed of astro- 
nomical silicates ([Weingartner fc Drainel[200lh . The grain 
mass density is fixed to 3.5g/cm'^. 

The dust optical properties: extinction and scattering 
opacities (Fig[T]), scattering phase functions, and Mueller 
matrices are calculated using the Mie theory. The resulting 
midplane optical depth in I band (0.81 ^m), from the star 
to the observer, is ranging from 1.22 x 10"^ to 1.22 x 10^ 
when the disc mass varies from 3 x 10~* to 3 x 10~^ M©. 
For simplicity, in the following, we will label the different 
models r = 10^, r ^ 10"^, r = 10^ and r = 10^. 
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Fig. 2. Mie scattering phase function (first element of the 
Mueller matrix, full line) for the 1 /Ltm silicate grains at a 
wavelength of 1 /xm. The dashed line represents the Henyey- 
Greenstein phase function for the same asymmetry param- 
eter g = {cos 9) — 0.63. 



In the simplifying case of Mie scattering, the matrix 
becomes block-diagonal with only 4 non-zero elements: 
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where the individual elements Sij only depend on the scat- 
tering angle and not on the azimuthal angle. The first ele- 
ment Sii, also known as the "phase function" is plotted in 
Fig [5] for the wavelength of 1 fj,m used to calculate scattered 
light images and polarisation maps. 

A spherical grain of 1 /zm at a wavelength of 1 fim is 
in the middle of a resonance region, where constructive 
and destructive interference within the dust grain results in 
phase (and polarisation) functions with strong oscillations. 
These effects are not observed in the case of a grain size dis- 
tribution (where the oscillations corresponding to different 
grain sizes are averaged) or in the case of more naturally 
shaped particles, like aggregates. However, we choose these 
dust properties as they represent a better test case for the 
radiative transfer codes. The oscillations in the elements 
of the Mueller matrix must be seen in the synthetic maps 
allowing a more detailed comparison between codes. For 
comparison an Henyey-Greenstein phase function with the 
same asymmetry parameter is over-plotted in Fig [51 

All the calculations presented in section U] are done as- 
suming anisotropic scattering and use the previously pre- 
sented Mueller matrix. To compare the results of the Monte 
Carlo codes with those obtained from discrete ordinate 
codes, the temperature structures and SEDs are also calcu- 
lated in the isotropic case {i.e. Sn is constant), using the 
opacities calculated with the Mie theory (section (5]). 

3.3. Maps and SEDs 

SEDs, images and polarisation maps are calculated at 10 
inclinations equally spaced in cosine, i.e. for cos(i) = 0.05, 
0.15, . . . , 0.85 and 0.95. This corresponds to inclination 



angles ranging from 18.2 to 87.1° from pole-on. Scattered 
light images and polarisation maps are calculated at 1 /zm. 
The pixel scale is 25.61 mas. pixel^^ {i.e. 251 pixels for a 
physical size of 900 AU at a distance of 140 pc). This is 
roughly a factor 2 smaller than the pixel scale of the WFPC 
and ACS cameras on-board the Hubble Space Telescope. 



4. Results 

4.1. Temperature structures 

Figures [3] and 2] show the temperature distributions cal- 
culated by the different codes, as well as the difference 
between codes. Overall, the agreement is very good with 
difference almost always smaller than 10%. 

Figure [3] presents the temperature along the disc mid- 
plane for the T = 10"^ and t = 10^ cases. Very close to the 
inner edge, where the disc is directly heated by the star, 
the agreement between codes is excellent with maximum 
differences of the order of 1%. At large radii (> 100 AU), 
the disc becomes optically thin at optical wavelengths in 
the vertical direction, and the midplane is heated by the 
stellar light that is scattered in the surface layers of the 
discs. In these regions, the peak-to-peak differences between 
codes remain below 1.5 and 3% for r = 10^ and r = 10^ 
cases, respectively. This shows that all codes deal smilarly 
with the redistribution of energy by anisotropic scattering. 
This is confirmed by the vertical cut at a radius of 200 AU 
(Fig m right panel) , where the differences are of the order 
of 1 %, except at the turnover point from optically thin to 
optically thick (the place where the temperature suddenly 
drops) where differences reach 5%. 

Differences in the radial temperature profile are larger 
between the inner edge and lAU, i.e. in regions were the 
stellar radiation does not penetrate, even via scattering. In 
these regions, the heating mechanism is the dust re-emission 
by the upper layers of the disc, which represent the most 
difficult case for RT codes. Nevertheless, the peak-to-peak 
differences between codes remain limited, smaller than 4 % 
for the r = 10"^ case. For the r = 10^ cases, differences are 
most of the time smaller than 10%, except in a very small 
regions between r-jn + 10^^ AU and rjn + 10^^ AU where the 
maximum difference is 20%. Differences remain also very 
small in the vertical direction as shown in the left panel of 
Fig SI They are below 5 % from the midplane up to the disc 
surface, where they become smaller than 1 %. 

4.2. SEDs 

The emerging spectral energy distributions for the r = 10^ 
and r = 10^ models are presented in Fig [5] for different 
inclinations ranging from an almost face-on (i = 18.2°) to 
an almost edge-on disc {i — 87.1°). 

The shape of the SED is strongly dependent on the in- 
clination, moving from a stellar photosphere plus disc ex- 
cess for low inclinations to a double-bumped SED, typical 
of very close to edge-on systems, when the stellar photo- 
sphere is obscured by the disc. The disc being optically 
thick at short wavelengths, the visible and near-IR stellar 
light is blocked and the emission in this wavelength range 
is dominated by the stellar scattered light coming from the 
disc. At longer wavelengths (> 10 — 12 pm) the dust emis- 
sion dominates, resulting in a steep positive slope and a 
double-bumped SED. Not surprisingly, the transition be- 
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Fig. 3. Radial temperature profile in the disc midplane. The left panel corresponds to the t = 10^ case and the right 
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dot-dashed lines the results of MCMax. The results of ProDiMo, RADMC and RADICAL are presented in Fig [H] (isotropic 
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tween these two characteristic SEDs also depends on the 
disc opacity. For instance, for an inclination of 75.5°, the 
star is still seen directly in the r = 10^ case, whereas it is 
already strongly obscured in the r — 10^ case. 

It is interesting to note that for the most edge-on case, 
the flux in the optical is almost independent of the disc 
opacity. Similarly, in the most optically thick case, there is 
only small differences in the SEDs at 81.4 and 87.1°. Indeed, 
when the star is completely obscured, the flux is dominated 
by stellar light that has scattered on the upper layers of 
the outer disc. This scattered light is mainly a function 
of the dust properties and of the scattering geometry but 
not of the optical depth in the line of sight. In these cases, 
optical and near-infrared photometry cannot be used to get 
estimates of the extinction by dust of the central object. 

At low inclinations, the characteristic 9.8 /Ltm amorphous 
silicate feature is seen in emission. At higher inclinations, 
for example i = 75.5° for the r = 10^ case, the feature is 
now observed in absorption over the continuum emission of 
the disc. It should be noted that, at very high inclinations, 
although the dip is roughly centred on the silicate feature, 
it is not associated to it as demonstrated by the exceptional 
breadth of the feature. The disc is now optically thick in 
the silicate feature but also in the adjacent continuum, and 
the absorption feature from the silicates vanishes. 

At long wavelengths (> 500 /xm), the disc become op- 
tically thin in most of its parts and the emerging flux no 
longer depends on the system inclination. 

Overall, the differences between codes are small in all 
cases, which is illustrated by the thinness of the grey en- 
velopes around the lines in Fig[5l representing the range of 



the results obtained by the different codes. Fig [6] gives a 
more detailed view of the differences between the codes for 
the various inclinations and optical depths. 

The left column show the results for the results for the 
T = lO'^ case. When the star is seen directly (i = 18.5 
and 75.5°), the agreement between the three codes is excel- 
lent, with peak to peak difference smaller than 5 % over the 
whole wavelength range. Closer to edge-on, differences re- 
main smaller than 10%, except at very short wavelengths 
(< 1 /im) where the contribution from scattered light is 
dominating the SED. As the wavelength becomes shorter, 
the scattering becomes more forward throwing and calcula- 
tions are very sensitive to the angular sampling of the scat- 
tering phase function of the codes. As a result, the agree- 
ment between codes becomes worse at short wavelengths. 
Differences remain smaller than 15% down to 0.2 /im and 
signiflcant differences are only observed around 0.1 /im. 

For the i = 81.4° case, the disc is seen at a grazing 
incidence and the optical depth from the star to the ob- 
server is varying strongly across the stellar disc, from w 2 
at the top of the stellar surface to w 200 at the bottom. 
In this case, using a point source for the star does not pro- 
vide the correct result, and special care must be taken to 
resolve the stellar photosphere. In the case of a uniformly 
radiating sphere, as presented in this benchmark, the ori- 
gin of the photon packets is uniformly distributed on the 
stellar surface, and a uniform distribution in the cosine of 
the angle between the photon direction and the normal to 
the surface at the point of origin is used to set the initial 
propagation direction of packets. 
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Fig. 5. Spectral energy distributions for i=18.2, 75.5, 81.4 and 87.1° (full, dashed, dotted and dot-dashed lines respec- 
tively). The left panel corresponds to the r = 10^ case and the right panel to the t = 10^ case. The lines show the 
average of the results of MCFOST, MCMax and TORUS and the grey envelope around each line represents the range of 
results obtained by the different codes. 



For the t — 10^ case (right panel of Fig[S]), the results 
are similar to the ones for the lower optical depth case, but 
with larger differences in the near- and mid-infrared. This 
part of the SED is one of the most challenging wavelength 
range for RT codes, where different contributions (thermal 
emission from the inner disc seen directly or through the 
outer disc, direct or scattered stellar light, scattered ther- 
mal emission from the disc) can dominate the emerging 
flux depending on the system geometry and dust opacities. 
Furthermore, most of the flux in this wavelength regime 



is coming from the inner edge of the disc and the output 
spectrum is very sensitive to the grid resolution adopted 
by the different codes. MCMax and MCFOST present very 
good agreement over all the wavelength range, including in 
the near- and mid-infrared regime, and for all inclinations. 
TORUS shows slightly larger differences, probably due to 
a lower spatial resolution at the inner edge. We note that 
the spatial resolution of the TORUS code is limited by the 
maximum cell depth in the AMR grid, which is currently 
set to 30, correponding to a dynamical range of 2'^'' (a limit 
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dictated by numerical precision in the photon path inte- 
grator). These differences are maximum around 10 /im and 
vary from 15 % in the low inclination case up to 30 % when 
the inclination is increasing. 

4.3. Scattered light images and polarisation maps 

Fig [7] presents the scattered light images of the most opti- 
cally thick disc for i = 69.5° and i — 87.1°. The synthetic 
maps clearly display the effects of the anisotropy of the 
scattering. The oscillations in the phase function (Fig ^ 
are directly observed in the maps. 

The three panels on the right present the flux obtained 
by the various codes, and the corresponding differences, 
along horizontal and vertical cuts in the images. The codes 
agree to within 10% where the flux is significant. TORUS 
shows slightly larger deviations that are due to a larger 
Monte carlo noise in the images. All codes predict the same 
oscillations as a function of the position, indicating that the 
implementation of the scattering phase function is correct 
in all codes. At greater radii, larger differences are observed 
due to the various grid geometries used by the different 
codes. For instance, some codes use a spherical grid whereas 
other codes use a cylindrical grid with a vertical cut-off in 
density. The sampling of the density structure at the outer 
edge of the model domain is then slightly different for each 
code. This results in systematic differences between codes, 
but only in the regions of the synthetic maps where the flux 
is extremely low. 

The vertical cut (#3) samples the disc "dark lane" , cor- 
responding to the optically thick midplane. In this region, 
the flux is dominated by photons that have scattered sev- 
eral times in the disc before reaching the observer. The very 
good agreement between codes in the dark lane indicates 
that all of them deal properly with multiple scattering. 

Fig [5] presents the polarisation maps for the same con- 
figurations as in the Fig. [7l The maps show complex pat- 
terns due to the strong variations of the elements of the 
Mueller matrix with the scattering angle. All codes pro- 
duce very similar maps (see horizontal and vertical cuts), 
at both inclination angles, with differences smaller than 5 
points of polarisation degree in the central regions of the 
maps, i.e. in regions where the flux is large enough to allow 
resolved polarisation measurement in actual observations. 
At greater radii, differences become slightly larger, but they 
remain limited to 10 points of polarisation degree. TORUS 
shows larger deviations, due to a larger Monte Carlo noise 
in the simulations, which biases the polarisation degree to- 
wards larger values. In regions where the polarised flux is 
large, the agreement between TORUS and the other codes 
is also very good. For the highest inclination, the results 
from TORUS are not shown due to a low signal-to- noise. 

5. Comparison with discrete ordinate codes 

All results presented thus far were obtained using Monte 
Carlo codes. To further compare the predictions of radiative 
transfer codes, we present, in this section, results obtained 
with discrete ordinate codes and discuss their limitations 
compared to those of Monte Carlo codes. Discrete ordinate 
codes solve the radiative equation along predetermined sets 
of directions. This integration can be performed with "long" 
or "short characteristics" and various schemes can be used 



to iterate between the temperature structure and specific 
intensity (or its moments), like the Accelerated Lambda 
Iteratio n (ALI) or Variable Edding ton Tensor (VET) meth- 
ods (see lMihalas fc Mihalaslll984 . 

The benchmark problem presented in this section is the 
same as in the previous sections but with the following mod- 
ifications: 

— the star is now considered as a point source, but with the 
same spectrum and luminosity as previously defined. 

— the scattering is assumed either to be isotropic (but with 
the same opacities as before) , either to be negligible (the 
scattering opacity is set to 0). 

These tests have been defined to allow a larger number of 
codes to reprodu ce th e calculations. They correspond to an 
extension of the |P04| benchmark to optical depths higher 
than 100, but without any further complexity. In this paper, 
they have been calculated by at least one of the previously 
tested Monte Carlo code and by two additional discrete 
ordinate codes and one additional Monte Carlo code. 

5.1. Code description 
5.1.1. ProDiMo 

ProDiMo is an acr onym for Protoplanetary Disk Model 
(jWoitke et al.l I2009D which consistently solves the chem- 
istry, the heating/cooling balance of the gas, the dust ra- 
diative transfer and the vertical stratification of protoplan- 
etary discs, mainly for the purpose of interpreting far IR to 
mm gas emission lines. 

For realistic gas models, it is essential to calculate the 
dust temperature structure in the disc as well as the trans- 
port of UV photons including scattering, which drive the 
photo-chemistry. Furthermore, radiative pumping by con- 
tinuum radiation changes the non-LTE population of atoms 
and molecules and have an important impact on the cooling 
rates. These strong physical couplings necessitate to solve 
a full 2D dust radiative transfer as one module in a global 
iterative procedure. It is this radiative transfer module in- 
side ProDiMo that participates in this benchmark test. Its 
basic task is to provide 2d(r, z) and J\{r, z) for the gas mod- 
elling - it is not meant for the interface to dust observations 
(SEDs, scattered light images etc.). 

ProDiMo solves the frequency-dependent 2D dust con- 
tinuum radiative transfer of irradiated discs by means of 
a simple, ray-based, long-characteristic method. From each 
grid point in the disc, a limited number of rays (here 172) 
are traced backwards, while solving the radiative transfer 
equation with isotropic scattering. The setup of the ray di- 
rections is critical for the optically thin parts of the disc, 
in particular at near IR wavelengths where the illumina- 
tion originates from small and far hot regions. This is done 
in a manual fashion in ProDiMo to ensure there are more 
rays pointing toward the hot inner regions than toward the 
cooler interstellar side. One central ray pointing toward the 
star is reserved and covers the solid angle occupied by the 
star. 

Instead of a treatment with a large number of wave- 
length grid points, ProDiMo uses a coarse wavelength grid 
{Xk\k^ 0, ...,K] (here = 24) from 100 nm to 1000 ^m, 
and treat the opacities, intensities and source functions 
with band means, e.g. Bk{T) = J^^ ^ B\{T) dX where 
AAfc ~ Xk — Xk-i- One radiation transfer iteration takes 
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Fig. 6. Differences in the SEDs obtained by the different codes. The disc opacity and inchnation are given on top of each 
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about 4 seconds for a low resolution 20 x 20 grid, and about 
70 seconds for a 70 x 70 grid on a single- processor 2.66 GHz 
Linux machine, which is comparable to the computational 
efforts taken to solve the disc chemistry. 

In order to solve the condition of radiative equilib- 
rium and the scattering problem, a simple A-type itera- 
tion is applied. The source functions are pre-calculated on 
the grid points and fixed during one iteration. During the 
ray-tracing, the opacities and source functions are interpo- 
lated from the grid point values. After having solved all 
rays from all points in all frequency bands, the mean in- 
tensities are updated and the dust temperatures are re- 
calculated. Without further accelerations, this A iteration 
converges only for problems up to a midplane optical depth 
of about 100. In ord er to a c celera te the convergence, we ap- 
ply the procedure of lAueil (|1984f ) known as "Ng" -iteration. 
This enables us to solve radiative transfer problems up to 
T = 10'^ . . . lO'', depending on the geometry of the model. 

For higher optical depths, we apply an approximate pro- 
cedure following an idea of CP. DuUemond which consists 
in reducing the dust density in the central midplane regions 
in the following way. For every vertical column (consider- 



ing the downward direction) we do not increase the dust 
density any further once a certain critical optical depth at 
1/im is reached (rcrit^lO), provided that the radial optical 
depth is also > Tcrit . With this "trick" , we can manage test 
problems up to r = 10^ with this code. 

5.1.2. RADMC 

RADMC is a Monte-Carlo based continuum radiative trans- 
fer code for 2-D axisymmetric configurations such as cir- 
c umstellar discs and e i ivelop es. The basic algorithm is that 
of lBiorkman fc WoodI ()2001[ ). but with a continuous deposi- 
tion of energy instead of the discrete deposition as described 
in the original paper. In this way the temperature profile is 
also smooth in the very optically thin regions of the model. 
The temperature corrections are not computed every time 
a photon package enters a cell and leaves some energy, but 
only if the energy deposited since the last temperature up- 
date is larger than some threshold value. In this way the 
not so cheap temperature update is done only when needed. 
Of course, the higher this threshold, the faster the code is 
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but the less reliable the result. Typically a threshold of a 
few percent is taken, meaning that if the energy deposited 
in the cell has increased by more than a few percent of the 
energy as it was during the last temperature update, then 
a new temperature update is done. Once the main Monte 
Carlo process is over, the spectra and images from RADMC 
are made using a post-processing step: the ray tracing pro- 
gram RAYTRACE uses the dust temperature and isotropic 
scattering source function computed by RADMC to calcu- 
late the formal solution of the transfer equation along rays 
through the model. This yields images at every discrete 
wavelength bin. By integrating over the images one obtains 
a flux at each wavelength, i.e. a spectrum. Care is taken 
to arrange the pixels of the images such that all the flux 



is captured, both from the very outer regions and from the 
very inner regions. This is done using "circ u lar irn ages" , de- 
scribed in detail in iDullemond fc Turollal (|2000[ ). Because 
the spectra and images are computed as a post-processing 
step, as opposed to the more classical photon collection dur- 
ing the Monte Carlo process itself, it is hard to include full 
non-isotropic scattering. To keep the flexibility to view the 
object from every angle and at every wavelength without 
having to repeat the RADMC run, one has to store the en- 
tire scattering source function S{r, 9, ii, (jj, A), where /i and (j) 
are the local directional coordinates. Such a 5-dimensional 
array is extremely large and requires of the order of a gi- 
gabyte of disc space, which is not very practical. Instead 
one could prescribe before calling RADMC at which angle 
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you wish to view the object, removing the need to store the 
source function also as a function of /i, cf), meaning we have 
a 3-D array which is much easier to store. This is done in 
several other codes in this paper. This is not implemented 
in RADMC. 

5.1.3. RADICAL-VET 

The code RADICAL is a classical discrete ordinate method 
for radiative transfer, i.e. it is not based on a Monte 
Carlo approach and is therefore completely deterministic. 
It is based on methods that are routinely used in mod- 
els of stellar atmospheres, with some adaptions. The al- 
gorithm is that of Variable Eddington Tensors (VET), 
which is a multi-dimensional version of the method of 



Variable Eddington Factors (VEF) described in the book by 
Mihalas & Mihalas (1984). A 1-D version of this method, 
with special application to the kind of continuum radia- 
tive transfer pro blems encountered in pro toplanetary discs, 
was described in lDullemond et all (|2002| ). For such 1-D ge- 
ometries the method is extremely accurate and efficient. It 
works well and converges quickly for optical depths rang- 
ing from small (<^Cl) to extremely large (;2> 10^). The 
RADICAL-VET code is a 2-D version of this algorithm. For 
2-D or 3-D geometries the method has some numerical diffi- 
culties related to the computation of the flux-mean opacity 
in regions of extremely low flux (e.g. the midplane of a pas- 
sive irradiated disc). In practice, however, these difficulties 
are not fatal, although they could lower the reliability of 
the method in such flow-flux regions. For further details on 
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Fig. 9. Radial temperature profile in the disc midplane (isotropic scattering). The left panel corresponds to the r — 10^ 
case and the right panel to the r — 10^ case. The full red lines represent the results of TORUS, the blue dashed lines, 
the results of MCFOST, the black dot-dashed lines the results of RADMC and the green dot-dot-dashed line the results 
of ProDiMo. Differences are relative to the average results of the Monte Carlo codes MCFOST, RADMC and TORUS. 



the V ET method the reader is referred to lDullemond et aTl 
()2002f ). even though this paper describes only the 1-D ver- 
sion of this method. 

5.2. Results and discussion 

The midplane temperature profiles are presented in Figs M 
and[TO]for the cases with isotropic scattering and no scatter- 
ing respectively. Convergence was not properly reached for 
some codes for the highest mass case with isotropic scatter- 
ing and we present the results for the r = 10^ case instead. 

The agreement between Monte Carlo codes is again very 
good in the isotropic case (Fig [9]) , with peak-to-peak dif- 
ference smaller than a few percents in the t = 10'^ case 
and smaller than 15 % in the r — 10^ case. Differences with 
discrete ordinate codes are significantly largers. 

This paper states an important verification test for 
the method of reducing the dust density implemented in 
ProDiMo and described in section 15.1.11 It shows that 
the error of this procedure in comparison to the results 
of the latest Monte Carlo codes combined with diffusion 
solvers is smaller than about 20 % in the midplane regions. 
We note however that the relative precision in general is 
much better than in the midplane. The average differencial 
between ProDiMo and MCFOST temperature structures, 
|TproDiMo - rMCFOSTl/rMCFOST is Smaller than 5 % (1 sigma- 
deviations) in all cases. 

The benchmark test has revealed three major problems 
quite typical for ray-based codes. First, the limited num- 
ber of fixed rays causes some small artifacts of the tem- 
perature structure in the optically thin distant midplane 



^ over 500 representative points in the grid 



- these problems can be solved by using a larger num- 
ber of rays which is, however, computationally expensive. 
Second, the precise temperature determination in the op- 
tically thick core of the disc is hard with simple discrete 
ordinate codes Hke ProDiMo. The numerical solution of the 
radiative transfer equation has always some discretisation 
errors superimposed, and in an optically thick situation it 
is the small difference I\ — B\ that determines the next 
temperature iteration. This eventually limits the quality of 
the "forecast" by the Ng-iteration, and disables the conver- 
gence for very optically thick problems. Third, the results 
close to the midplane suffer from the very large gradients 
present close to the inner boundary, and the results depend 
on the numerical details, e.g. how to interpolate the source 
function. 

Comparisons with the RADICAL code illustrates some 
of the difficulties of the VET method at very high optical 
depths. RADICAL overestimates the midplane temperature 
in the central regions of the disc by about 20 % and 40 % for 
the low and high mass cases respectively (Fig [TOl) . The 2-D 
geometry introduces a number of difficulties which make the 
2-D VET algorithm less stable and less reliable than its 1-D 
version. The main problem is the choice of discrete angu- 
lar coordinates at each grid point. For the formal transfer 
RADICAL-VET uses the method of Short Characteristics. 
The choice of the angular distribution of these short char- 
acteristics is essential for t he reliability of the result. In 
iDullemond fc TuroUal (|2000D a good choice was described, 
but in the end no choice is perfect and the reliability of the 
results may depend on this. Another difficulty is that the 
VET method uses flux-mean opacities which are the gener- 
alisations of Rosseland mean opacities. By using flux-mean 
opacities instead of Rosseland mean opacities we may ex- 
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Fig. 10. Radial temperature profile in the disc midplane (no scattering). The left panel corresponds to the lowest disc 
mass of 3 10^^ M0, and the right panel to the highest disc mass 3 10^^ M©. The corresponding optical depths at 1 //m are 
lower than in the previous cases because we fixed the scattering opacity to zero. The red full lines represent the results 
of ProDiMo, the blue dashed lines the results of MCFOST and the black dot-dashed lines the results of RADICAL. The 
Monte Carlo code MCFOST is taken as reference. 



pect to get the true result instead of just an approximate 
result. But in regions where the flux is extremely small, such 
as in the midplane of an extremely optically thick disc, tak- 
ing the average of the opacity based on a quantity that is 
nearly zero is dangerous and may lead to large errors. In 
spite of these caveats the VET algorithm gives reasonable 
results for most problems. 



6. Summary 

We have presented solutions for the continuum radiative 
transfer in high-opacity circumstellar discs. The problems 
have optical depths up to 10^ and include anisotropic 
scattering. They represent realistic configurations for discs 
around low mass stars and validate the use of the codes to 
model current and future observations of discs. 

We have compared the results of four independent 
Monte Carlo codes for the temperature structure, SEDs, 
scattered light images and/or polarisation maps, in the case 
of anisotropic scattering. Overall, the agreement between 
codes is very good. In the most optically thick case, SEDs 
agree within 20 % over almost all of the wavelength range. 
Differences become larger only at wavelengths shorter than 
0.2 /zm for edge-on configurations, i.e. when the fiux is ex- 
tremely low and not observed in practice. Pixel-to-pixel dif- 
ferences in high-resolution scattered light images remain 
limited to 10 % and the polarisation maps do not differ by 
more than 5 points of polarisation degree in regions where 
the polarisation can be effectively measured by observa- 
tions. Each observation (SED, image or polarisation map) 
was reproduced by at least three of the codes, providing 
robust solutions to test other RT codes. 



The benchmark problems represent challenging test 
cases for RT codes. The convergence of Monte Carlo meth- 
ods alone become extremely slow for the most optically 
thick cases and specific numerical schemes are required 
to efficiently compute the temperature structure, emerging 
SEDs and images, for instance combining a Monte Carlo 
approach with ray-tracing and/or diffusion approximation 
methods. 

Comparisons between Monte Carlo codes and discrete 
ordinate codes, in the cases with isotropic scattering or 
without scattering, show relatively large differences, that 
increase with optical depth. Ray-tracing and VET meth- 
ods were successfully compared to Monte Carlo methods 
up to moderate optical depths (ry = 100, lP04l ) but the 
convergence of such codes seem to become delicate in the 
test cases presented here. They provide good approximate 
solutions but must be used with care at high optical depths, 
when the goal is to perform detailed comparisons with ob- 
servations. 
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